#include <cusp/io/matrix_market.h>
#include <cusp/array2d.h>
#include <cusp/multiply.h>
#include <cusp/print.h>
#include <cusp/transpose.h>
// #include <cusp/qr.h>
//#include <cusp/qrhh.h>
// #include <cusp/qrhhgemv.h>
#include <cusp/qrhhb.h>
#include <cusp/detail/frobenius_norm.h>
#include <cusp/detail/profiler.h>

#include <iostream>
#include "../../performance/timer.h"

int main(void)
{
    // create a simple example:  It is crucial to use column-major storage for cublas calls
    cusp::array2d<float, cusp::host_memory, cusp::column_major> Ahost, Qhost, Qthost, Rhost, Rthost, AoutHost;
    cusp::array2d<float, cusp::device_memory, cusp::column_major> A, Q, R;

    cusp::io::read_matrix_market_file(Ahost, "/project/csstaff/inputs/Matrices/MatrixMarket/RECTANGULAR/well1033.mtx");

    double m = (double)Ahost.num_rows;
    double n = (double)Ahost.num_cols;
    double mflops = 4.0 * n * ( m * (m - n) + n*n/3.0) / 1000000.0;
    // print A
    std::cout << "Matrix rows " << Ahost.num_rows << " cols " << Ahost.num_cols << " pitch " << Ahost.pitch << "\n";
//    cusp::print(Ahost);
    A = Ahost;  // copy to device

    std::cout << "Perfoming QR decomposition (Givens)..." << std::endl;
//    cusp::qr(Ahost,Qhost,Rhost);  

//    Qhost = Q; Rhost = R;
//    cusp::multiply(Qhost, Rhost, AoutHost);
//    cusp::print(AoutHost);

    std::cout << "Perfoming QR decomposition (Householder)..." << std::endl;
//    cusp::qrhh(A, Q, R);
//    cusp::qrhhgemv(A, Q, R);
    cusp::qrhhb(A, Q, R, 20);   // Warm up GPU
    {
    timer t;
//    cusp::qrhh(A, Q, R);
//    cusp::qrhhgemv(A, Q, R);
    cusp::qrhhb(A, Q, R, 20);   // Time GPU
    float timeD = t.seconds_elapsed();
    std::cout << "Device time: " << timeD << " s " << mflops/((double)timeD) << " mflop/s " << std::endl;
    }

//    Qhost = Q; Rhost = R;   // copy back to host
    // cusp::transpose(Qhost,Qthost);
//    cusp::multiply(Qhost, Rhost, AoutHost);
//    cusp::print(Ahost);
//    cusp::print(AoutHost);
    CUSP_PROFILE_DUMP();
    
    return 0;
}

